This Page Is Inserted by IFW Operations 
and is not a part of the Official Record 



BEST AVAILABLE IMAGES 

Defective images within this document are accurate representations of 
the original documents submitted by the applicant. 

Defects in the images may include (but are not limited to): 

• BLACK BORDERS 

• TEXT CUT OFF AT TOP, BOTTOM OR SIDES 

• FADED TEXT 

• ILLEGIBLE TEXT 

• SKEWED/SLANTED IMAGES 

• COLORED PHOTOS 

• BLACK OR VERY BLACK AND WHITE DARK PHOTOS 

• GRAY SCALE DOCUMENTS 



IMAGES ARE BEST AVAILABLE COPY. 



As rescanning documents will not correct images, 
please do not report the images to the 
Image Problem Mailbox. 



WORLD INTELLECTUAL PROPERTY ORGANIZATION 
International Bureau 




PCX 

INTERNATIONAL APPUCATION PUBUSHED UNDER THE PATENT COOPERATION TREATY (PCT) 



;51) International Patent Classification : 
GOIV y48 



Al 



(11) International Publication Number: 
(43) IntemaUonal Publication Date: 



WO 00/60380 

12 October 2000 (12.10.00) 



[21) IntemaUonal Application Number: PCr/IBOO/00353 

(22) IntemaUonal Filing Date: 27 March 2000 (27.03,00) 



I April 1999 (01.04.99) 



GB 



(71) Applicant (for GB oniy): SCHLUMBERGER LIMITED 
[NLOJS]; 277 Park Avenue, New York. NY 10172 (US). 

(71) AppUcant (for AT AU DE DK ES ID IE IT KR NO NZ RU 
only): SCHLUMBERGER TECHNOLOGY B.V. [NUNL]; 
Parfcstraal 83-89, NL-2514 JG TTie Hague CNL). 

(71) AppUcant (for BR MX only): SCHLUMBERGER SURENCX) 
S A. [PA/PA]; 8 Calle Aquilino dc la Guaidia. Pananw City 
(PA). 

(71) Applicant (for BF BJ CF CG CI CM GA GN GWML MR NE SN 
TD TO only): PETROLEUM RESEARC:H AND DEVEL- 
OPMENT N.V. [NL/NL]; De Ruyterkade 62, Willemstad, 

Curacao (AN). 

71) Applicant (for CA only): SCHLUMBERGER CANADA LIM- 
ITED [CA/CAl; 24th floor. Monenco Place, 801 6fli Avenue. 
S.W., Calgary. Alberta T2P 3W2 (CA). 



(71) Applicant (for CN only): SCHLUMBERGER OVERSEAS 
S.A. PPA/PA]; 8 Callc Aquilino dc la Guardia, Panama City 
(PA). 

(71) Applicant (for FR only): SERVICES PETROUERS 
SCHLUMBERGER P^l/FR]; 42. me Saint Dominique, 
F-75007 Paris (FR). 

(71) Applicant (for all designated States except US): SCHLUM- 

BERGER HOLDINGS LIMITED (-/-]; P.O. Box 71. Craig- 
muir Clhambcrs, Road Town, Tortola (VG). 

(72) Inventors; and 

(75) Inventors/Applicants (for US oniy): BRIE, Alain [FR/JPJ; 
4-2-12 Dcnenchofu Ota-ku, Tokyo 145-0071 (JP). 
ENDO. Takeshi (JP/JP]; Grace CHiikazawa 603. 6-19-15 
Sagamihara, Sagamihara-«hi, Kanagawa-ken 229-0031 
(JP). VALERO, Henri-Pierre [FR/JP]; Urvan Kakishima II 
#502. 6-19-9 Haramachida, Machida-shi, Tokyo 194-0013 
(JP). PISTRE, Vivian [FR/US]; 2323 Long Reach Drive 
#10303, Sugar Land. TX 77478 (US), SAIKI, Yoshiyuki 
[JP/JPl; Cities Shine No. 10. 202, 25-16 Susukino-cho. 
Sagamihara-shi. Kanagawa-ken 229-1113 (JP). 

[74) Agent: HYDEN, Martin, Douglas; Schlumberger-Riboud Prod- 
uct Ontre, Dept. Brevets, 26, fue de la Cavee. B-92142 
Oamart (FR). 



(81) Designated States: AE, AL. AM, AT, AU. AZ. BA, BB, BG. 
BR. BY. CA, CH, CN, CR. CU. CZ, DE. DK. DM, EE. 
ES. n, GB, GD, GE, GH. GM, HR, HU. ID, IL, IN. IS, JP. 
KE, KG, KP, KR, KZ, LC LK. LR. LS, LT. LU. LV, MA, 
MD, MG, MK, MN, MW. MX, NO, NZ, PL, PT, RO, RU, 
SD, SE, SG, SI. SK, SL. TJ, TM. TR. TT, TZ. UA, UG, 
US, UZ, VN, YU, ZA, ZW. ARIPO patent (GH, GM. KE. 
LS, MW, SD, SL. SZ, TZ, UG, ZW), Eurasian patent (AM. 
AZ, BY, KG, KZ. MD. RU, TJ, TM), European patent (AT. 
BE, C:H, CY, DE. DiC ES, FI. FR, GB, GR, IE, IT, LU. 
MC, NL, PT, SE), OAPI patent (BF, BJ. CF. CG, a, CM, 
GA, GN. GW. ML. MR. NE, SN. TD. TG). 



Published 

With international search report. 



(54) TiUc: PROCESSING SONIC WAVEFORM MEASUREMENTS 



200 
A 100 



sot 



305 3)0 319 320 325 330 335 3^ 345 3S0 3SS 



(57) Abstract 



A method of determining the sonic slowness of an underground formation firom sonic measurements, includes the steps oft (i) 
obtaining sonic waveforms in the underground formation; (ii) determining, from the sonic wavefomis, multiple values of at least one 
parameter related to the sonic slowness of the formation together widi an estimate of the error in each value; and (ui) usmg the estimate of 
the error in each value to select a parameter value related to slowness as representative of the formation. The sonic waveforms are typically 
obtained by logging an interval of a borehole which runs through the fonnation with a tool which output sonic waveform measiuements. 
The tool can be nin on wireline or coiled tubing or die like, or alternatively can be a logging while drilling tool located m a drill stnng 
being used to drill the borehole. 
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PROCESSING SONIC WAVEFORM MEASUREMENTS 
TECHNICAL FIELD 

5 The present invention relates to methods for processing sonic waveform measurements, 
particularly sonic waveform measurements made for the purpose of characterising properties of 
underground formations. The invention in particular relates to methods for determining the best 
value for a parameter that has been determined in a number of different ways. 

10 BACKGROUND ART 

It has been know for some time that it is possible to determine properties of underground 
formations using measurements of acoustic/sonic waves that have passed through the formations. 
The basic technique comprises placing a tool comprising a spaced sonic source and receiver in 

15 the borehole and using the source to generate sonic waves which pass through the formation 
around the borehole and are detected at the receiver. Sonic waves can travel through rock 
formations in essentially two forms: body waves and surface waves. There are two types of body 
waves that travel in rock: compressional and shear. Compressional waves, or P-waves, are waves 
of compression and expansion and are created when rock formation is sharply compressed. With 

20 compressional waves, small particle vibrations occur in the same direction the wave is travelling. 
Shear waves, or S-waves, are waves of shearing action as would occur when a body is struck 
from the side. In this case, rock particle motion is perpendicular to the direction of wave 
propagation. The surface waves are found in a borehole environment as complicated borehole- 
guided waves which come from reflections of the source waves reverberating in the borehole. 

25 The most common form or borehole-guided, surface wave is the Stoneley wave. Figure I shows 
a series of sonic waveforms such as would be recorded in a borehole from a monopole 
(omnidirectional) source with the first arrivals of the compressional (P), shear (S) and Stoneley 
(St) waves at the receiver marked. In situations where dipole (directional) sources and receivers 
are used, an additional shear/flexural wave propagates along the borehole and is caused by the 

30 flexing action of the borehole in response to the dipole signal from the source. The flexural wave 
typically travels at about the same speed as the shear wave, slower than the compressional wave. 
(It is to be noted that sonic waves will Also travel through the fluid in the borehole and along the 
tool itself. With no interaction with the formation, these waves carry no useful infomiation and 
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can interfere with the waveforms of interest if they have similar propagation speeds.) A further 
phenomenon affecting sonic measurements in boreholes is that fluid waves, Stoneley and 
flexural are dispersive; that is their velocity is a function of frequency. 

The speeds at which these waves travel through the rock are controlled by rock mechanical 
properties such as density and elastic dynamic constants, and other formation properties such as 
amount and type of fluid present in the rock, the makeup of the rock grains and the degree of 
intergrain cementation. Thus by measuring the speed of sonic wave propagation in a borehole, it 
is possible to characterise the surrounding formations by parameters relating these properties. 
The speed or velocity of a sonic wave is often expressed in terms of 1/velocity and is called 
"slowness". Since the tools used to make sonic measurements in boreholes are of fixed length, 
the difference in time (AT) taken for a sonic wave to travel between two points on the tool is 
directly related to the speed/slowness of the wave in the formation. 

15 An example of a tool for use in a borehole for sonic measurements is the DSI tool of 
Schlumberger which is shown schematically in Figure 2. The DSI tool comprises a transmitter 
section 10 having a pair of (upper and lower) dipole sources 12 arranged orthogonally in the 
radial plane and a monopole source 14. A sonic isolation joint 16 connects the transmitter 
section 10 to a receiver section 18 which contains an array of eight spaced receiver stations, each 

20 containing two hydrophone pairs, one oriented in line with one of the dipole sources, the other 
with the orthogonal source. An electronics cartridge 20 is connected at the top of the receiver 
section 18 and allows communication between the tool and a control unit 22 located at the 
surface via an electric cable 24. With such a tool it is possible to make both monopole and dipole 
measurements. The DSI tool has several data acquisition operating modes, any of which may be 

25 combined to acquire (digitised) waveforms. The modes are: upper and lower dipole modes 
(UDP. LDP) - waveforms recorded from receiver pairs aligned with the respective dipole source 
used to generate the signal; crossed dipole mode - waveforms recorded from each receiver pair 
for firings of the in-line and crossed dipole source; Stoneley mode - monopole waveforms from 
low frequency firing of the monopole source; P and S mode (P&S) - monopole waveforms from 

30 high frequency firing of the monpole transmitter; and first mouon mode - monopole threshold 
crossing data from high frequency firing of the monopole source. 
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One way to determine compressionaU shear and Stoneley slownesses from these measurements 
is to use slowness-time-coherence (STC) processing. STC processing is a full waveform analysis 
technique which aims to find all propagating waves in the composite waveform. The processing 
adopts a semblance algorithm to detect arrivals that are coherent across the array of receivers 

5 and estimates their slowness. The basic algorithm advances a fixed-length lime window across 
the waveforms in small, overlapping steps through a range of potential arrival times. For each 
time position, the window position is moved out linearly in time, across the array of receiver 
waveforms, beginning with a moveout corresponding to the fastest wave expected and stepping 
to the slowest wave expected. For each moveout, a coherence function is computed to measure 

10 the similarity of the waves within the window. When the window time and the moveout 
correspond to the arrival time and slowness of a particular component, the waveforms within the 
window are almost identical, yielding a high value of coherence. In this way, the set of 
waveforms from the array is examined over a range of possible arrival times and slownesses for 
wave components. STC processing produces coherence (semblance) contour plots in the 

15 slowness/arrival time plane. Regions of large coherence correspond to particular arrivals in the 
waveforms. The slowness and arrival time at each coherence peak are compared with the 
propagation characteristics expected of the arrivals being sought and the ones that best agree 
with these characteristics are retained. Classifying the arrivals in this manner produces a 
continuous log of slowness versus depth. For dispersive waves, the STC processing is modified 

20 to take into account the effect of frequency. As the output of STC processing is a coherence plot, 
the coherence of each arrival can be used as a quality indicator, higher values implying greater 
measurement repeatability. When processing dipole waveforms, one of the coherence peak will 
correspond to the flexural mode but witii a slowness that is always greater (slower) than the true 
shear slowness. A precomputed correction is used to remove this bias. 

25 

To compensate for variations in measurements due to the borehole rather than due to the 
formation a series of measurements are made across an interval in which the formation 
properties are expected to vary littie, if at all. In its simplest form, the interval corresponds to the 
extent of the receiver array, and the waveforms at each receiver station measured for a given 
30 firing of a source ("receiver array" or "receiver mode" or "Rec"). In simple STC processing, all 
receiver stations are considered. In multishot STC processing (MSTC), sub-arrays of receiver 
stations witiiin the receiver array are considered, for example a sub-array of five receiver stations 
in a receiver array of eight receiver stations (other numbers or receiver stations in the sub-array 

3 
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can be used depending on requirements). In this case, the same formation interval corresponding 
to ihe extent of a five receiver station sub-array can be measured several times as the tool is 
logged through the borehole, the five stations making up the sub-array being selected at each 
source firing to measure the same formation interval. Another approach, known as '"transmitter 

5 mode" or "pseudo-transmitter array" ('Tra/') takes waveforms from sequential source firings as 
the transmitter passes along the interval to be measured. In order to compensate for the 
movement of the tool between measurements, an effectively stationary receiver station or sub- 
array must be used. This can be achieved by changing the receiver station considered so that its 
position in the borehole is effectively stationary as the transmitter is moved through the interval. 

10 Borehole compensation ("BHC") can be achieved for P and S mode results by processing 
receiver array and pseudo-transmitter array waveforms and averaging the results. 

Thus it will be appreciated that, with the different acquisition modes of a tool such as the DSI, 
and the different processing modes available, it is possible to obtain multiple determinations of a 

15 slowness or AT in a given interval of borehole. For example, it is possible to acquire waveforms 
for shear slowness determination using two dipole modes, (upper dipole and lower dipole). and 
one monopole mode (P and S mode), and to process each measurement in receiver mode, 
transmitter mode and borehole compensated form resulting in a potential nine separate 
determinations of shear slowness for that interval, each of which can give a different result. The 

20 problem is therefore to determine which slowness estimate can be considered to give the best 
indication of the shear slowness of the formation in that interval. 

DISCLOSURE OF INVENTION 

25 The present invention provides a method of determining the sonic slowness of an underground 
formation from sonic measurements, comprising: (i) obtaining sonic waveforms in the 
underground formation; (ii) determining, from the sonic waveforms, multiple values of at least 
one parameter related to the sonic slowness of the formation together with an estimate of the 
error in each value; and (iii) using the estimate of the error in each value to select a parameter 

30 value related to slowness as representative of the formation. 

The sonic waveforms are preferably obtained by logging an interval of a borehole which runs 

through the formation with a tool which outputs sonic waveform measurements. The tool can be 

4 
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run on wireline or coiled tubing or the like, or alternatively can be a logging while drilling tool 
located in a drill string being used to drill the borehole. 

The multiple values of the parameter related to sonic slowness can include mutliple 
5 determinations of monopole and/or dipole compressional and shear (including flexural) 
slowness* and Stoneley slowness for the formation. Where the tool used to obtain the waveforms 
comprise an array tool, the multiple determinations can include transmitter and receiver mode 
measurements and borehole compensated measurments. The processing technique used is 
preferably a slowness-time-coherence technique. 

10 

The processing technique can include amongst its inputs, zoning information derived from the 
waveform measurements and indicating general features of the formation type being measured. 
The zoning information can be obtained from a basic compressional slowness estimation, 
typically based on digital first arrival determination, and can include broad formation slowness 
IS classifications such as fast, slow, very slow and extremely slow. Such broad classifications can 
be based on predetermined cross plots of the ratio of compressional and shear slownesses against 
measured slowness for know lithologies. Other zoning information can be the presence of 
closely spaced bed boundaries (thin beds). 

20 Other inputs to the processing technique can include parameters relating to the borehole or well, 
such as hole diameter, mud type and predetermined formation features. 

The specific processing technique applied can vary according to the zoning or parameter inputs. 
For example, where zoning information shows relatively thick beds, full array STC-type 

25 processing (including dispersive processing) can be applied; if zoning information shows thin 
beds, high resolution, sub-array multishot STC-type processing can be applied. The processing 
also preferably includes tracking of slowness measurements along the interval so as to allow 
individual slowness measurements to be associated with changes in particular components of the 
waveform (P-waves, S-waves, etc.). The tracking can also make use of the zoning information 

30 indicating where major changes occur, and delineating homogenous beds, each with associated 
semblance error bars. 



5 
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The output of the processing step is a series of slowness estimates and the tracking can also 
include the estimation of error in the value of slowness. An estimate of error can be provided by 
the statistical semblance processing of the slowness measurements. The total error is a 
combination of the semblance error and the uncertainty in tracking determination. The final step 
selects the slowness with the minimum error (possibly modified by other predetermined 
selection rules) and outputs this as the slowness for the interval. For example, the mean error for 
the interval can be used as the basis for selection. Also, a level by level determination of the 
variation of the actual error from the mean within the interval can be provided as a further 
output 

BRIEF DESCRIPTION OF DRAWINGS 

Figures 1 shows a plot of sonic waveforms in time indicating arrivals of various components; 
Figure 2 shows a prior art sonic logging tool; 

Figure 3 shows a flow diagram outiining the main features of a method according to one 

embodiment of the invention; 

Figure 4 shows a simple compressional slowness log; 

Figure 5 shows a VpA^s vs compressional slowness cross plot; and 

Figure 6 shows a slowness log with multiple slowness determinations. 

BEST MODE FOR CARRYING OUT THE INVENTION 

Referring now to Figure 3, the method summarised therein takes as its inputs tiie raw waveforms 
100 obtained from a borehole logging tool such as die DSI tool from Schlumberger, zoning 
information 110 and processing parameters 120. The raw wavefonns 100 are essentially digital 
signals comprising the receiver station output in time for a given acquisition mode of the tool 
(cf. Figure 1). To obtain the zoning information 110. the raw waveforms 100 are pre-processed 
using a digital first arrival detection (DEAD) method such as is described in more detail in 
International Patent Application No. W097/28464 to derive a crude compressional slowness 
(DTCO) log vs depth 102 an example of which is shown in Figure 4. The DTCO log 
information is analysed 104 to determine the general zones present in the interval of interest. 
This is achieved by applying predetermined thresholds to the DTCO log and squaring the output 
so as to broadly classify the log into zones of fast, slow, very slow and extremely slow formation 

6 



wo 00/60380 PCT/IBOO/00353 

with sharp transitions between the zones (line DTCOsq in Figure 4), and is called '^macro- 
zoning". Where the zones are changing frequently over relatively short distances (thin beds, for 
example < 2ft thickness), the log is also indicated as having "micro-zoning". The macro- and 
micro-zoning information is output to a AT computation process 130 (micro-zoning information 
5 can be used to indicate bed boundaries and to select multishot STC processing which is 
described in more detail below) and (for macro-zoning information)to a parameter selection 
process 112. 



The parameter selection process 112 determines from physical well parameters 114 a series of 
10 processing parameters (see examples given in the table below) and provides these as a 
processing parameter input 120 to the AT computation process 130. 



Well Parameter 




Tool Tvne 




Formation Type 


Select from automatic using zoning input or 
manual selection based on VpVs-DTCO 
crossnloi overlaid with formation tvoe regions. 


Hole Diameter 


Bit size, casine ID or caliper 


Borehole Status 


Ooen Hole or Cased Hole 




Processing Parameter 




Dispersion Curve 


Slowness dispersion curve for the mode to be 
evaluated. _ 


Integration Time Window 


Window size for coherence calculation, few 
cvcles of lowest freauency component 


Slowness Lower Limit 
Slowness Upper Limit 
Slowness Step Parameter 


Stan, stop and increment along slowness axis 
of S/T plane for STC processing. 


Time Lower Limit 
Time Upper Limit 
Time Step Parameter 


Start, stop and increment along time axis of 
S/T plane for STC processing. 


Search Band Width 
Search Band Offset 


Define limits of S/T plane to be processed. 


Slowness Width 
Time Width 


Define mask for peak search. 


Semblance Threshold 


Minimum semblance value for peaks. 


Filter Length, Upper Limit, Lower 
Limit 


Processing frequency band. 


' Mud Global Parameters 


Mud type: Brine, WBM, OBM, Slow OBM; or 
calculate from mud type (brine, emulsion, 
invert emulsion) salinitv, density, temperature. 



7 
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It will be appreciated that not all of these parameters might be needed, or that default values 
might be acceptable in some cases. The particular parameters, how they are derived and their 
significance will depend on the particular processing scheme used. There are other parameters 
which are optional such as optimal frequency selection for processing such as is described in 
5 more detaU in US 5,587,966. 

The parameter "Formation Type" can have two main methods of selection, both of which rely on 
the zoning step 104. In the first, the DTCO log data for a given level is plotted on a scale which 
is banded into predefmed formation types. In Figure 5, DTCO data from two depths are 
indicated by + and x. For the purpose of illustration, these are cross plotted with VpA^s and the 
predefined formation types: fast (F), intermediate (I), slow (S). very slow (VS) and extremely 
slow (XS) superimposed on the plot The predefined formation types are derived from other log 
data obtained in a number of locations. The data from depth + fall in a fast formation type while 
those from depth x fall in a slow formation type. Thus different parameters will be selected for 
processing the data from these two depths. It will be appreciated that, normally, there will be no 
Vp/Vs dau available and it is only position on the DTCO axis which is used to decide formation 
type. An alternative method is to output DTCOsq directly as the formation type. 

The AT computation process 130 comprises two main elements, an integrated slowness 
20 determination process 132 and a tracking process 134. The integrated slowness determination 
process 132 provides a Slowness/Time Coherence (STC) methods for analysing the waveforms. 
These methods can include one or more of STC, multishot STC (MSTC), fast STC (FSTC) and 
dispersive STC ODSTC). These different methods are summarised below. In the preferred 
example, the major processing method utilised is DSTC. 

25 

Slowness Time Coherence (STC) Processing 

STC processing employs full waveform analysis to find all propagating waves in an array of 
sonic waveforms. The processing computes a peak vector that identifies all of the arrivals in the 
waveform data from a sonic array tool such as the DSL The peak vector consists of seven 
30 elements associated with a peak coherence value in the slowness/time plane. Two main steps are 
carried out in the STC technique: 

8 
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1. The coherence o(T,S) is calculated for all reasonable values of Time and 
Slowness in the S/T plane. o(ZS) serves as a measure for whether the data 
includes an arrival at time T and slowness S. 

2. Then, the surface a(T,S) is searched for local maxima with a peak finding 
process. 

The coherence measure is calculated by : 



Where Ec(T,S) is the coherent energy of the normalised waveforms (square of the summed 
10 waveform values), Ei(T,S) is the total energy of the normalised waveforms (sum of the squared 
wavefonn values), and M is the number of waveforms, i.e. receivers (receiver stations). 

Ec(T,S) and Ei(T,S) are calculated over a specified time window, usually selected to span several 
cycles of the main fi^uency component expected in the waveforms. In Ec(ZS), the wavefonns 
15 are stacked in time assuming the arrivals move out linearly across the receiver array. For any 
number of waveforms, the coherence measure a(T,S) will be between 0 and 1. Values 
approaching 1 indicate the presence of a highly coherent arrival at time Tand slowness 5. 

The coherence a(ZS) is calculated within a restricted search band of the S/T plane defined by 
20 the slowness lower and upper limits and time lower and upper limits. For a given slowness S 
between the lower and upper slowness limits, the time (7) in the search band must obey the 
following rule: 



where ZO is the T/R spacing of the tool. Toff is the time offset of the search band, and is the 
25 width along time of the search band. 

This rule describes a band of times moving di^igonally across the S/T plane from the lower, left- 
hand comer to the upper, right-hand comer. For Toff equal to 0, this band is centred on the 
slowness-time line (r = 5 • ZO) . 




S • Z0 + % - nand /2 < r < 5 • Z0+ % +7^„^ /2 



30 



SUBSTITUTE SHEET (RULE 26y. 
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The point (Tpk, Spfc) is defined as the peak coherence in the S/T plane search ask if it satisfies the 
following: 

1 • t5r(Tpk, Spk) > 0.25 (Threshold value) 

A rule of thumb is that values greater than 2/M indicate the presence of an arrival. 
If 8 wavefomis are being input, the value of 0.25 can be assumed. 
2. a (Tpk, Spk) is a local maximum within the given peak (specified by the slowness 

width an time width): 
<^(Tpk^Spk)> p{T„,.Sj for all T„ and S„ 

such that 

r.W/2<7^<7-hW/2 

where T^ask is the time width of the peak mask, and Smask is the slowness width of 
the peak mask. 

IS The following peak matrix of seven elements are output to seven arrays at every depth 
processed: 
At peak coherence: 





1. 


o(Tpb Spt) 


coherence value 




2. 


Spk 


slowness 


20 


3. 




time 




4. 


Ec(Tph Spk) 


coherent energy in decibels 




At peak energy: 






5. 


0(T^ Spk) 


coherence 




6. 


Te 


time 


25 


7. 


Ec(T^ Spk) 


coherent energy in decibels 



The slowness is the same as for peak coherence. 



Signal-to-noise ratio is estimated using the following algorithm: 

1. The peak-to-peak amplitude in the signal window is determined at each receiver 
30 and the amplitudes are averaged to get an estimate of the signal. 

2. The peak-to-peak amplitude in the quite window is determined at each receiver 
and the amplitudes are averaged to get an estimate of the noise. 

10 



5 
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3. The signal-to-noise ratio estimate is then simply calculated using the following 
equation: 5 / = 20 • Log 1 ^{signal I noise) . 

There are two main depth reference points in the STC processing (or DSTC processing, see 
5 below), waveform depth and computation depth. Waveform depth is the depth assigned to all 
waveforms in the array at a particular depth level and is assumed to be the mid point between the 
fired transmitter and the first receiver in the receiver array. All waveforms in the array are 
assigned this depth. The computation depth is the depth assigned to the S/T plane and the peak 
matrix computed by one STC computation. When a receiver array mode is selected, this depth is 
10 the raid-point of the receiver array used in the computation. Thus if all eight receivers are used, 
the depth is the mid-point between receivers 1 and 8. If a sub-array of receivers 5 - 8 is used, the 
depth is the mid-point between receivers 5 and 8. When a transmitter array mode is used, the 
depth is the mid-point of the transmitter positions forming the pseudo-array. 

15 Dispersive STC (DSTC) Processing 

Certain modes of sonic waves are dispersive, for example the dipole flexural shear mode. This 
means that different frequencies propagate at different phase slownesses, waveshapes change as 
they propagate across the array, and measurement of shear slowness requires a dispersion model. 
Processing the flexural mode with standard, non-dispersive processing such as STC (described 

20 above) does not give foirnation shear slowness without the use of model-based corrections. 
Dispersive STC (DSTC) processes the flexural mode dispersively so no correction for dispersion 
is required. In STC, waveforms in the window are backpropagated to the reference receiver 
along lines of constant slowness. This constant slowness is the same slowness as the moved-out 
window and the calculation can be done entirely in the time domain. The semblance of the 

23 backpropagated waveforms is calculated, high semblance values indicating that the 
backpropagated waveforms are nearly alike which, in turn, indicates that the backpropagating 
slowness is correct. DSTC performs backpropagation in a different way. The waveforms are 
backpropagated by a flexural mode dispersion curve equal to that of the data. The dispersion 
curve is pre-computed from the flexural mode model and for speed the backpropagation is 

30 performed in the frequency domain. The shear slowness of the best-fitting curve is relumed by 
DSTC. DSCrr also outputs an error bar for the uncertainty of the slowness determination due to 
noise. Further details of DSTC can be found in US 5,278,805, 
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DSTC processing has the ability to be applied even to non-dispersive waves such as monopole 
compressional or shear.. Since the first step required for DSTC processing is the calculation or 
selection of an appropriate dispersion curve, all that is required is a dispersion curve that 
represents a non-dispersive wave, Le. a flat "curve". Consequently, DSTC is the preferred 
processing scheme for the present invention and will be indicated as such in the following 
description. However, it will be appreciated that in many cases, STC processing could equally be 
applied. 

Multishol STC (MSTC) Processing 

MSTC attempts to improve the resolution of the DSTC full array calculation by processing 
results from sub-arrays and oversampling the interval of interest. For example, in an eight 
receiver tool, it is possible to improve resolution by forming smaller sub-arrays of, for example, 
five receivers. While this reduces the amount of data available from each sub-array, the fact that 
the tool has multiple opportunities to sample the region with different sub-arrays of the same 
size means that this can be compensated in the redundancy of the data. MSTC does this using the 
following procedure: 

1. All of the sub-arrays which span the interval of interest are identified. 

2. Each sub-array is processed using DSTC (an S/T plane is generated for each sub- 
array). 

3. The S/T planes for the sub-arrays are shifted in time to account for the difference 
in transmitter-receiver spacing. 

4. The semblances from the various sub-arrays are averaged at every value of time 
and slowness to provide and average S/T plane. 

5. The average S/T plane is searched for slowness peaks which are returned in a 
peak matrix. 

The peak matrix output by MSTC contains three elements: 

- the slowness of the coherence peak, 

- the coherence of the coherence peak, and 

- the time of the coherence peak. 

The slowness, time and coherence are simply the slowness, time and coherence of each peak 
found in the average S/T plane. 
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Selection of MSTC as the processing is made when the zoning information 110 indicates thin 
beds (typically <2ft thickness). A single-shot DSTC compulation by MSTC is the same as an 
STC computation made by DSTC. Further details of MSTC can be found in US 4,809,236. The 
general processing methodology preferred for the present invention is described in further detail 
5 below. 



Fast STC (FSTC) Processing 

Fast STC (FSTC) processing was developed for use with LWD sonic measurement techniques 
but is potentially an alternative to the standard STC processing described above. In FSCT the 
10 raw waveforms are filtered to give a complex signal type with both real and imaginary parts. 
Because of the increase in data arising from this, the data are coarsely resampled to reduce the 
amount of data and only the real part of the complex signal data are saved for display purposes. 
The basic processing other than the filtering and sampling are as for STC. Further details of 
FSTC can be found in US 5.594,706. 

15 

Whichever processing method is used, the output is typically more than one slowness for each 
depth, e.g. corapressional slowness, shear slowness and Stoneley slowness. Figure 6 shows a log 
with more than one slowness (e.g. compresional, shear, Stoneley) indicated for each level. It is 
generally desired to monitor the development of these properties along the well and so it is 

20 necessary to decide which slowness is which from depth to depth in the well. This can be 
difficult in cases where the slownesses are similar in value or cross each other, or are missing 
because of bad data at a given depth. Commonly used techniques for connecting peaks in a 
continuous curve are called labelling and some use a logical decision tree, level by level. More 
sophisticated techniques use a weighted cost function over a short interval to select an optimum 

25 path. The particularly preferred approach is to use a two-step process 134, first joining the peaks 
that correspond to the same waveform arrival ("u-ack search"), and second identifying the u-ack 
with a name ("classification"). 



To perform the mick search, the sequence of measurements that are sufficiently continuous are 
30 enumerated and their probabilities computed. In order to know how to associate the different 
peaks with the different track built in this way, all possible hypotheses are computed and the 
corresponding probabilities derived. With this information, and by applying Bayes theorem, the 
a posteori probability can be computed and the hypothesis with the most probable value selected 

13 
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to build a track. Figure 6 shows point on the log connected to form tracks. Classification is 
achieved by fonning different hypotheses based on a propagation model to give all possible 
classifications of tracks built, and associating a priori probabilities with these hypotheses. 
Knowing the data and the a priori probability model, the a posteori probabilities of the different 
5 hypotheses can be computed and the most probable hypothesis considered to classify the tracks 
136. Since the classification is statistical in nature, the error in the determinations can be 
estimated and ou^ut as a quality indicator in each case 138. The zoning information 110 is also 
input to this determination since it helps identify where transitions between formation types 
might cause sudden changes in the slowness curves which otherwise might be misinterpreted. 

10 

The track search methodology is as follows: Where is the set of all measurements up through 
the depth ik and Z(k) is the measurement at the same depth, if at depth k n peaks are recorded, the 
measurement Z(k) is composed of n slownesses, n times and n semblances. 2^ consists of all 
measurements up to and including the considered depth: 

15 Z^={Z(0),Z(l),...,Z(*)} 

If g/* is the i-th global hypothesis for the measurements 2*, up through depth t, each g/" gives 
one possible complete explanation of the measurements up through the considered depth. Thus 
the sequences of measuiements continuing over depth (tracks) are completely specified and the 
measurements not taken into account are considered as false alarms. This hypothesis can be 

20 decomposed as two parts qi(k) and qi(i)(k). The first hypothesis, qi(kK specifies at depdi k, how 
measurement Z(k) is assigned to tracks and how tracks at depth k-l are linked with tracks at 
depth Jk. The hypothesis, qi(i)(k). specifies the l(i) hypothesis for measurement for measurements 
through to depth Z*"^. This decomposition can be rewritten as follows: 

Z^' =Z^"^ uZ{k) 

25 qi^qi(i)^M 

This shows that it is possible recursively to generate all hypotheses for all measurements up to 
the current depth, k. The tracking algorithm works by forming hypotheses about measurements 
and evaluating their probabilities using prior probabilities model and Baye's rule. The 
probability of the global hypothesis 9/ based on all available measurements Z^ can be expressed 
30 as: 

14 
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fz(t:j,f.z*-'XJ(*Kv^'"'K4i 



This equation is the general equation of the tracking algorithm without any simplification 
assumptions. The denominator Q is a normalising constant obtained by summing the numerator 
of the equation over /. Conceptually, the tracking search part is completely defined by this 
equation. Each time a new frame of data Z(k) is acquired, the hypotiiesis, qi(k), at this depth is 
generated and appended to the hypothesis for the measurements up through the previous 
depth, qi(i)^'^ and finally the prior probability and data likelihood models are used to compute 



10 The last term of the numerator in the equation, P(qt^\2^'^), represents the probability of the 
parent global hypothesis which is considered to be known from the previous depth. The two 
other terms of the numerator equation can be obtained according to Kurien, T., 1992, Issues in 
the design of practical multi-target tracking algorithms,, Multitarget-multisensor tracking: 
applications and approaches,, Artech House. The term P(qi(k)\qi(i^'\2l''^ ) is defined as 

p[gr\Oi.Ni.Z^''y(h,^^^^^ 

where 0,% the modelorder, gives the number of tracks at one depth k implied by qf(k) and Oia) is 
the modelorder implied by ^ Ni is the sum of all the tracks observed at the depUi k and the 
number of false alarm measurements implied by qi(k), gAOi,Ni) is considered as the r-th way to 
classify Ni measurements as false alanns and arrivals up to 0/ arrivals. In fact the current 

20 modelorder O, and the number of measurements Ni determine the number of possibilities in 
which Ni measurements are classified as false alarms and arrivals up to a maximum of O, 
arrivals. The probabihty P(gr\OiMi) is computed with the assumption that each (9, arrival is 
independent of any others and of the false alarms. Supposing that for one frame of O, arrivals 
each has a probability p to appear at the current depth and a probability l-p not to appear. The 

25 false alarms are modelled as independent and their number is modelled as Poisson distribution 
with mean 1. Considering K of the 0, arrivals present in the classification gr(Oi,Ni) then tfiere are 
K arrivals and Ni -K false alarms and the probability is: 

15 
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hs(0i,0i(i)) is a function that gives the s-th possible way to assign 0/ tracks at depth k with the Oia) 
tracks in the previous depth k-1. This function gives an enumeration of the different ways to link 
tracks in the current depth with those in the previous depth. Considering that all the combinations 
are equally probable, for y/(Oi,Oi(i)) possibilities, then 

The classification of measurements as arrivals or false alarms, gAOiyNO^ and track assignment, 
hs(Oi,Oi(i)), with the probability models P(hs\Oi,Oi(i)) and P(gAOi,Ni)\Oi,Ni) give a complete 
determination of P{qi(k)\qi(i^'\2!''^ ), 

The last term of the concept equation is computed as follows: Given data Zj(k) at depth k 
decomposed into a continuous part C(k) and a discrete part N(k): Z(ifc) = {c(A:),A^(ik)}. 

Considering C(k) and N(k) are independent, P^^ikjjfif ^Z^^^ ^ can be rewritten as: 
P{z{k\qlz^-' y P[c{k\qf .Z^-' y[N{k\qlz^-'^ ) 
15 where p\N{k]^qf .Z*"^ j= 1 if Ni=N(k) and 0 otherwise. 

q/^ completely describes the sequences of measurements which belong to the same tracks. To 
know if the sequences of measurements enumerated as "sufficiently continuous" to qualify as 
real tracks or not, the probability likelihood model for the continuous part of the data 

20 P^C{kj^qf , Z*"^ j must be provided. 



10 



C(k) is formed by the slowness-time co-ordinate of each peak modelled as a 2D measurement 
vector. Assuming that the prior probability model considered is that the sequences of slowness 
and time are Gaussian random process, the correlation of slowness and time along depth is 
25 modelled as the output of an ARMA filter with known coefficients. To compute this probability, 
Kalman filter theory is used. Scharf, L, 1991, Statistical signal processing: Wiley, is used to 
implement the ARMA filter via a Kalman state representation and Harvey, A., 1994, Time series 
models: MIT Press^ to implement the Kalman update and prediction equations. Given the prior 

16 
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continuity model for a track /, and using successive conditioning, Kalman filter equations allow 
the computation of PfC^I/j which is the probabihty of continuous measurements from depth 0 to 

k: 



5 Each of the terms in the equation can be computed using Kalman prediction and update 
equations. At each depth Kalman prediction gives the 2D Gaussian probability distribution of 



the track is computed by evaluating the Gaussian distribution at the slowness-time point 
specified by C(i). Knowing this information, the 2D Gaussian probability distribution of 
10 slowness-time at depth i+l is predicted using C(i) and die result of updating of the state and 



needed to complete the likelihood computation above. 

By applying this methodology to the measurement data, measurements can be linked together on 
15 a depth-by-depth basis to form tracks. Thos measurements that cannot be so linked are rejected. 
It is then necessary to assign each track to a physical attribute, for example monopole 
compressional arrivals, monopole shear arrivals, etc. The approach used in this invention is to 
first make certain assumptions, for example that the first arrival track will be monpole 
compressional followed by monpole shear followed by Stonely. Given these assumptions, the 
20 data making up the tracks is then tested to determine the probability that a track fits these 
assumptions and are classified according to which assumption shows the highest probability. The 
general methodology for this determination is as follows: 

is the track data at depth k and contains the following information: called the modelorder, 
25 that gives the number of tracks identified at the k-ih depth, the triplet (P, /, r), the slowness, time 
and semblance of each track, and some linking information that provides how a track at depth A- 
is linked with the detected tracks at depths kA and /:+! . h/ is thej-th hypothesis at depth k. The 
index ; corresponds to the number of hypotheses in the particular case, for example in the case of 
monopole the value is 3 (for the three modes: compressional, shear and Stoneley). The posterior 
30 probability /// given the data at depth A-1 and all previous depths is defined by the following 
equation where (.) denotes {D^'\ D^): 




the slowness-time given the past measurements 



The likelihood of the data C(i) belonging to 



covariance matrix via the Kalman update equation. This is the final term 
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■ /fo'|H).(.)V(wjK.)) 

*!(.)) 

This considers that the data likelihood for the data at the current depth depends on the previous 
data only through the hypothesis at the current depth, which means .(.))= ). 

These two equations are the main equations of the classification. 

5 

The hypotheses for a given depth depend upon the number of tracks (m) at this depth. To assign 
the probabihties some hypotheses at the current depth will be taken as consistent 

with the hypotheses at other depths and other will not. For example, if the classification 
classifies one track at depth kA as a component A and classifies a track at depth as the same 

10 arrival A, this gives the most probable chance to classify the detected track as arrival A 
compared to any other. Therefore, P(H^\Ht'^) ^K \f the hypotheses Hj and H^'^ are consistent 
over the depth and zero otherwise. The value K is determined as follows: for M consistent 
hypotheses found, K=l/M. Thus K is linked with the number of consistent hypotheses found. In 
this case, all the hypotheses are taken into account as equally probable and it is assumed that the 

15 track search part has made no mistakes when linking measurements over the depth. Given a 
finite probability of error in the track search part, PiHjhnt'^ ) = p when the hypotheses are 
consistent over depth and l-p otherwise, with 0 < 1-/? < 1. In this case it is necessary to set the 
value of p inside the algorithm. 

20 The compulation of the hkelihood P(D^\Hf) assumes that the tracks are independent of each 
other, the model used to compute PiD^Hj) considers that each track is a slowness-time point in 
the 2D slowness-time plane which is Gaussian distributed. The prior probability of this 
slowness-time point depends upon the classification assigned to that point by Hj, If the number 
of tracks equals one, then the P(D^'\Hj')is a Gaussian distribution with mean and variance given 

25 by /// and evaluated at the slowness-time point given by D^'. For more than one arrival, the 
assumption is that the independency of tracks and P(D^\h/) is the product of all Gaussian 
distributions evaluated. 

18 
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It will be appreciated that the tracking and classification methodologies described above can be 
applied to sonic logging data independently of the other processing steps described here. 

The slowness logs 136 together with the appropriate error bars 138 are output to a finalisation 
process 140 which computes and selects the best compressional, slowness, the best shear 
slowness and the bets Stoneley slowness out of those available 142, e.g.: 
L Compressional Slowness: 

P&S Rec., Tra., or BHC 

UDP Rec.» Tra., or BHC 

LDP Rec.» Tra., or BHC 

2. Shear Slowness: 

P&S Rec, Tra., or BHC 
UDP Rec. Tra., or BHC 
LDP Rec, Thl, or BHC 

3. Stoneley slowness: 
MST Rec., Tra-, or BHC 

The process selects the log with the minimum error bar. 

The best slowness is computed from one or more datasets using the following steps: 

1. A mean error is computed over the whole interval over depth for every candidate. 

2. The mean error is computed for each processing mode. i.e. Rec. or Tra. If the 
candidate is processed by BHC mode, the mean error values from both processing 
modes are averaged. 

3. Absent value resident in the error bar logs is regarded as a large error in the error 
computation. 

4. A candidate is selected which has the minimum mean error. 

5. Slowness is computed level by level based on the error bars of the selected 
candidate. 

6. If an error at a level is bigger than a threshold, absent is assigned as a slowness at 
the level (no-log condition). If the candidate is processed by BHC mode, the 
procedure looks over both errors of receiver and transmitter mode. If their errors 
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are similar, the mean slowness is assigned as the slowness, otherwise the one with 
the smallest error is chosen. 
7. Final error bar logs associated with the selected slowness logs are derived in the 
same way 144. 

5 

The general multi-shot processing methodology preferred for the present invention is as follows: 



For the case of multiple shots at depth ^: 
Ai is the waveform taken into account to process MSTC. 
10 I is the number of a sub-array (/ = 7, AO- 

Fi is the ST plane computed for the sub-array i. 



A frame is represented as a function of two parameters: F/ = F, (r,5) where t is the time and S the 
slowness. For N sub-airays, iV = 1, 5 for example, five frames F = {Fj, F2, F3, F4, F5} are 
15 obtained where F is the family of frames at depth Z. 

The object of the multiple shot is to compute 

P'stack = X ^ 

But as there is a shift in time between two consecutive frames, it is necessary to take into 
20 account this time shift before computing the MSTC processing. 



Between frame F, and f ,.7 ,the tool moves a known distance d = Az.. For a frame F,(r,5). to shift 
this array before computation of the multiple shot, it is necessary to apply the following 
operation: for the point >i(//, 5/) located in the ST plane, it is necessary to shift it by: 
25 tti =ti-d'S{ti). 

For the totaUty of the frame, after time shift, this gives: 

Fi{t.S)^Fi{t-{i-l)d-S{tls). 

This equation, called the time shift equation, explains how to make a shift in time between one 
frame and the succeeding frame. Conceptually, to stack to ST plane of each array it is necessary 
30 to apply this correction. 

20 
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Starting with the equation above, it is possible to rewrite it as a general equation: 

Fi\t.S)^ Fiit.S)^ S{t - (/ - iV • S{t\s) 
where * is the convolution operator on the variable /. 

5 To resample the frame to stack it with others, the shah function is used, defined as: 

For a signal s(t) it is possible to rewrite this signal as: 

Sj{t)^s{t)shuh 

Thus the signal s(t) is discretised with a sampling rate r 

10 

Rewriting the equation of a frame shifted in time and after resampling gives: 



FirM=F;{t,S)shah\ 




Developing this equation gives: 

Fi^rM = F;{t.S)shah 

= [Fi (r, 5) * S{t - (/ - l)d S{t))]shah 

15 This defines the transformation to be applied on each frame before stacking. 

Knowing the transformation that it is necessary to apply on a frame i\ it is possible to write the 
general equation of multi-shot processing. For a given level Z, this is: 

Fstack =Fdt,S)+f^ p4F,{t,S)*S{t^{i-l)dS{t)Mi-nT). 

20 This equation is the general multi-shot equation independent of the medium (homogenous or 
heterogeneous). To avoid the problem of increase computation time due to resampling in time 
and time shift, an alternative is proposed that is conceptually similar to the multi-shot equation. 
The idea is to include the time shift problem during the processing of MSTC. Presently, when 
computing MSTC, one receiver is selected as reference for all of die processing. In the present 

21 
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case, the reference receiver of each sub-array is changed before computing the semblance of this 
sub-array. The preferred approach is to use as a reference the first receiver of the first sub-array, 
the second receiver of the second sub-array and so on. Each time the semblance for each of these 
sub-airays is computed. This takes into account the time shift between the different sub-arrays, 
5 so it is possible to stack the semblances immediately after computation, the time shift between 
the different sub-arrays being taken into account. 

After the time shift correction has been made, and the semblance for each sub-airay computed, it 
is necessary to stack the ST planes by using an algebraic mean to compute the final ST plane at 
10 one considered depth. This is achieved by computing: 

f=l 

where Ff is the ST plane computed for one sub-array, and N is the number of sub-arrays used to 
compute the final multi-shot result 

15 Efficient processing of multi-shot data requires that the following be determined: 

1 . The number of sub-arrays used to compute the multi-shot. 

2. Selection of good sub-arrays for computing multi-shot. 

3. Automatic detection of the reference receiver for each sub-array. 

20 For a tool with R receivers, Rsb is the number of receivers in a sub-array. This number is selected 
by the user and is linked to the resolution to be obtained. By definition, its lower limit is 3. The 
number does not change the actual processing, just the physical interpretation of the result of the 
processing. A^,^, is the maximum number of sub-arrays used to compute the multi-shot. This 
value is equal to the number of measurement depths through which the tool is moved to make 

25 the multi-shot measurements and can be defined as: 

In order to know which receivers it is necessary to use to compute the multi-shot, a *'geometrix 
matrix" is implemented. When computing a multi-shot, it is necessary to change the reference 
30 receiver for each sub-array when computing the ST plane to take into account the time shift 

correction. For Nsh sub-arrays of Rgb receivers, the geometrix matrix is defined as: 

22 
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The first receiver of the first sub-array is equal to the value of sub-array necessary to compute 
the multi-shot. Therefore, r« = Nsb^ the second receiver of this sub-array is r^+i and so on. Each 
column of the matrix can be computed as: 

Where / is the number of the receiver inside the sub-array. The reference receiver is always the 
same as the first receiver of the first sub-array. Thus if this value is computed at the beginning of 
the calculation, all of the parameters for the STC processing are fixed for this computation and it 
is not necessary to redefine the reference receiver for each sub-airay. Furthermore, the geometrix 

10 matrix allows easy identification at each frame the receivers that are activated and deactivated at 
that depth. Consequently, it is merely necessary to multiply the receiver status by the current 
vector to know which receivers are to be used, thus simplifying data management. All that is 
needed is to know whether or not a receiver is active for a given frame, since data from all 
receivers is calculated at the beginning and the results used by circular permutation of the 

15 column of the matrix for each new firame. 

It will be appreciated that this approach can be used for any multi-shot processing and not just 
in accordance with the present invention. 

20 

INDUSTRIAL APPLICABLmV 

The present invention finds application in the field of characterising underground formations 
surrounding borehole such as in the oil and gas industry. 
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CLAIMS 

1 A method of determining the sonic slowness of an underground formation from sonic 
3 measurements, comprising: 

(i) obtaining sonic waveforms in the underground formation; 

(ii) determining, from the sonic waveforms, multiple values of at least one parameter 
related to the sonic slowness of the formation together with an estimate of the error in 
each value; and 

10 (iii) using the estimate of the error in each value to select a parameter value related to 

slowness as representative of the formation. 

2 A method as claimed in claim 1, wherein the multiple values of the parameter related to 
slowness include multiple values of compressional, shear and/or Stoneley slownesses. 

15 

3 A method as claimed in claim 1 or 2, wherein the sonic waveforms are obtained using a 
tool having one or more transmitters and an array of receivers 

4 A method as claimed in claim 3, wherein sonic waveforms are obtained in transmitter 
20 mode, receiver mode and/or borehole compensated mode. 

5 A method as claimed in claim 3 or 4, wherein the sonic waveforms are derived from 
monopole or dipole measurements. 

25 6 A method as claimed in any preceding claim, wherein the step of determining the 
parameter related to sonic slowness includes the seep of determining a general formation 
type parameter from the sonic waveforms. 

7 A method as claimed in claim 6, wherein the general formation type parameter is derived 
30 from a compressional slowness log that has been squared to fit into predetermined bands. 

8 A method as claimed in any preceding claim, wherein the step of determining the 
parameter related to slowness further includes selecting processing parameters relating to 
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well properties, formation properties and waveform processing parameters, and applying 
these parameters to the processing of the sonic waveforms. 

9 A method as claimed in claim 8, wherein the processing parameters comprise a first set of 
S parameters that are selected by a user, and a second set of parameters that are selected or 

calculated on the basis of the first set of parameters. 

10 A method as claimed in claim 9, wherein the number of parameters in the first set is 
relatively small compared to the number of parameters in the second seL 

10 

11 A method as claimed in any preceding claim, wherein the waveforms are processed using a 
slowness time coherence technique. 

12 A method as claimed in claim 11, wherein a multishot technique is used in formations with 
15 thin beds 

13 A method as claimed in any preceding claim comprising linking the values of the 
parameter related to slowness along the interval of interest into tracks which describe the 
development of a particular slowness parameter in that interval. 

20 

14 A method as claimed in claim 13, further comprising assigning identiHers to each track 
which identify the slowness mode to which the track data relates. 

15 A method as claimed in claim 14, wherein the estimate of error in each value is determined 
25 for each track, 

16 A method as claimed in claim 15, wherein, for a given slowness mode with multiple 
candidates, the candidate with the minimum mean error is selected as representative of that 
slowness mode. 

30. 

17 A method as claimed in claim 16, further comprising outputting error logs for the selected 
candidate logs. 
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